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NUMERICAL ANALYSIS OF THE TRANSIENT RESPONSE OF 


ABLATING AXISYMMETRIC BODIES INCLUDING 
THE EFFECTS OF SHAPE CHANGE 

By Stephen S. Tompkins, James N. Moss, Claud M. Pittman, 
and Lona M. Howser 
Langley Research Center 

SUMMARY 

The differential equations governing the transient response of an ablating axisym- 
metric, orthotropic body have been derived for fixed points in a moving coordinate sys- 
tem. These equations have been expanded into finite-difference form and programed for 
numerical solution, with an implicit technique, on a digital computer. Numerical results 
compare favorably with exact solutions. 

Several applications of the analysis are discussed. These applications demonstrate 
the significance of a salient feature of the analysis, that is, the ability to analyze the 
effects of changes in body geometry. This feature was used to obtain satisfactory agree- 
ment between numerical and experimental results for an ablating teflon sphere and a 
small test specimen exposed to a high-intensity laser beam. Although the analysis is 
primarily for a single-layer material, a multilayer material can be successfully approx- 
imated under certain conditions by a single-layer system. 

INTRODUCTION 

One-dimensional ablation analyses have been used extensively to study the thermal 
response of heat shields subjected to aerodynamic heating. However, for heat shields 
with large curvature or large heating-rate variations over the surface, the assumptions 
of one-dimensional heat flow no longer apply, and an accurate description of the thermal 
response requires an ablation analysis for multidimensional heat transfer that includes 
the effects of changes in heat-shield geometry. 

References 1 to 6 are examples of two-dimensional thermal analyses presently 
available. The outstanding features of the analyses of references 1 and 2 are that both 
consider axisymmetric bodies of anisotropic materials and both use stable, implicit 
numerical methods to solve the heat-conduction equation. However, neither analysis 
considers mass transfer. Mass transfer is included in the analysis of reference 3 in 


addition to the indicated features of references 1 and 2. The analysis of reference 3 is 
formulated for fixed nodal points in a fixed coordinate system and often requires inter- 
polation at the boimdaries because the nodal points and the boundaries do not always 
coincide. This interpolation can lead to inaccuracies. The analysis of reference 4 is 
similar to that of reference 3 except that it uses a moving coordinate system (which 
eliminates interpolation) and a conditionally stable, time-consuming explicit formulation 
as opposed to the stable, time-saving implicit formulation used in references 1 to 3. 

Only references 5 and 6 consider the effects of shape change on the thermal response of 
the heat shield. However, reference 5 does not consider an anisotropic material, and 
both references 5 and 6 use explicit methods to solve the governing equations. 

Collectively, references 1 to 6 consider many of the significant physical character- 
istics of an axisymmetric ablating heat shield and also demonstrate the desirable method 
of solution, that is, an implicit method. However, no single reference incorporates all 
these features into one analysis. 

This paper presents and discusses a transient two-dimensional ablation analysis 
which incorporates all the significant characteristics of the physical problem considered 
collectively in references 1 to 6. The analysis has the following features: (1) the abla- 
tion material is considered to be orthotropic with temperature-dependent thermal prop- 
erties; (2) the thermal response of the entire body is considered simultaneously; (3) the 
heat -transfer and pressure distribution over the body are adjusted to the new geometry 
as ablation occurs; (4) the governing equations and several boundary- condition options 
are formulated in terms of generalized orthogonal coordinates for fixed points in a 
moving coordinate system; and (5) the finite -difference equations are derived and solved 
implicitly. 

The accuracy of the analysis presented in this paper is demonstrated by compari- 
sons between numerical results and exact solutions for simplified conduction problems. 
Selected examples and test data are shown to demonstrate the utility of the analysis and 
the importance of the coupling between body-shape change and the heating-rate and pres- 
sure distributions. 

The finite -difference equations have been programed for solution on a high-speed 
digital computer. The program has a plotting routine which can display the shape of the 
ablating body at any time during the calculation. 

SYMBOLS 


A =-^— defined by equation (10) 

Xb 

Ac constant in oxidation equation corresponding to specific reaction rate 
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Cj,Dj J 


coefficients in equation (28) 


Ag constant in sublimation equation 

a constant to adjust equations (C24), (C25), and (C26) to correct form for 

Cartesian coordinates 


B = hjhgk 
Be 


Bs 

C 

Cp 

H 

AHc 

AHs 


hi,h2,h3 


i 

K 

k 

L 


M 



constant in exponential of oxidation equation corresponding to activation 
energy 

constant in exponential of sublimation equation 

oxygen concentration by mass 

specific heat 

total enthalpy 

heat of combustion 

heat of sublimation 

coordinate scale factors (eqs. (2)) 

order of reaction (eq. (11)) 

reaction-rate constant for oxidation (eq. (15)) 

thermal conductivity 

number of stations in x-direction 

molecular weight of gas 

molecular weight of oxygen 



m,n 

m 

“O2 

mg 

P 

P\v 

‘ic 

,net 

‘Inet 

R 

^cyl 

^stag 

r 

S 

T 

Tb 

t 


integers 
mass loss rate 

rate at which oxygen diffuses to surface 
mass loss due to sublimation 

exponent of pressure in sublimation equation (eq. (17)) 
wall pressure 

convective heating rate to nonablating cold wall 

hot -wall convective heating rate corrected for transpiration (eq. (13)) 

net heating rate to surface including combustion, sublimation, and surface 
reradiation (eq. (19)) 

radiant heating rate 

radius of curvature of base curve 

cylindrical radius from axis of symmetry to base curve 
stagnation-point radius of curvature 

exponent of radius in sublimation equation (17); spherical coordinate 

number of stations in y-direction 

temperature 

temperature of body to which back surface radiates 
temperature at finite-difference station (m,n) 
thickness of heat sink 
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free-stream velocity 


Cartesian coordinates (see fig. 2) 
curvilinear coordinates (see fig. 1) 
length of base curve 
absorptance 

weighting factor for transpiration effectiveness of mass loss due to 
combustion 

weighting factor for transpiration effectiveness of mass loss due to 
sublimation 

either 0 or 1 depending on whether transpiration or ablation theory is used 

material thickness 

emittance 

angle between R and Rcyi (fig- 1); spherical coordinate 
mass of char removed per unit mass of oxygen 
dimensionless curvilinear coordinates, equations (4) 
density of material 
Stefan-Boltzmann constant 
time 

angle of rotation about axis of symmetry, figure 3 

angle between axis of symmetry and normal to surface, figure 1 



Subscripts: 


c combustion 

e edge of boundary layer 

m,n integers 

L last station in x-direction 

max maximum 

o original, value at previous time 

S last station in y-direction 

stag stagnation-point condition 

w wall condition 

x,y coordinates 

dimensionless coordinates 
Superscripts: 

' condition along x = L 

" condition along y = 0 


ANALYSIS 
Physical Model 

The analysis considers an axisymmetric ablating body exposed to aerodynamic 
heating; this body is composed of a single orthotropic material of varying thickness with 
temperature-dependent thermal properties. (See fig. 1.) Although the analysis con- 
siders a single-layer material, the analysis of a multilayer material, such as a charring 
ablator, can be successfully approximated with the present analysis under certain condi- 
tions. (See ref. 7, for example.) 


6 



Surface normal 



Two coordinate systems are used to study the thermal and ablative response of the 
heat shield. One is a curvilinear coordinate system with x,y coordinates (fig. 1), which 
is used to determine internal temperature distributions. A stationary base curve located 
at the back surface of the ablator establishes the x-axis. 

The second coordinate system (fig. 2) is used to define the heat-shield exterior 
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geometry which changes with time as a result of ablation. This coordinate system, with 
w,z coordinates, is a Cartesian system with the origin fixed at the original stagnation 
point on the heat shield. All the geometric parameters needed to compute changes in the 
stagnation heating rates and the heating-rate and pressure distributions over the surface 
are defined in this system. The equations that quantitatively relate the internal temper- 
atures with the changes in the heating and pressure over the surface of the heat shield 
are given in the following sections. 

Groverning Differential Equations 

The differential arc length ds in the curvilinear coordinate system in figure 3 is 
(ds)2 = hj(dx)2 + h2(dy)2 + h3(d(;o)2 (1) 

where the scale factors are 


hi 



(2a) 


h2 = 1 (2b) 

h0 = Rgyi + y cos 0 (2c) 


The curvilinear coordinate system should conveniently describe any axisymmetric 
body geometry of interest. However, if the Cartesian coordinate system is required, 
care must be taken to use unity scale factors. Unity scale factors are obtained by 
assigning the values of 1 to Rcyi, °° to R, and ir/2 to 9 . (See eqs. (2).) 


For an axisymmetric body, the governing time-dependent heat-conduction equation 
with variable coefficients is (in fixed coordinates) 


1 

8 Ai2h3 

^1^2^13 

8x1 hj 


8T\ . 8 [ hih3 
8^ ho 


8x 


+ 



8T 


(3) 


If equation (3) were expressed in finite- difference form, it would describe the tem- 
perature variation at fixed stations in a fixed coordinate system. To maintain a fixed 
number of stations in a layer which changes thickness with time, it is necessary to 
change the locations of the stations and to interpolate to determine the temperatures at 
the new location after each time step in the calculation. This procedure not only 
increases the time required to perform the computations but also introduces a small 
error in each step of the calculation. An alternative to this procedure is to transform 
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Figixre J.- Basic cixrvllinear coordinate system. 

the equation to a coordinate system in which the stations remain fixed and the coordinates 
themselves move to accommodate changes in the surface location. 

This transformation can be made by introducing a moving coordinate system which 
is defined by the following relations: 

^ ^ and 77 = ^ (4) 

Ab o 

In this system, the outer surface remains fixed at 77 = 1 and all other stations remain 
at fixed values of 77. 

Before equation (3) can be transformed to the ^,77 coordinates, derivatives with 
respect to x and y in terms of ^ and 77 must be determined. These derivatives 
are given in the following equations: 
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( 5 ) 


D - J-— 1^ — 

"ix ~ ^ ^ 0X 8?7 ~ x^j 8 ^ ' 6 3x 877 


and 


^A = I_L 

8y 8y 8| 8y 877 6 877 


( 6 ) 


Because 77 is a function of 6 which is also a function of time, the time derivative on 
the right side of equation (3) becomes 


IH— A. 

dr ~ dr 8| dr 877 dr 


(7a) 


The derivatives of | and 77 are, from equations (4), 




= = 0 
8t Xjj 8t 

77 86 

dr 86 dr 6 8 t 




(7b) 


A change in 6 is given by 


Therefore, 


6-6 


£6 ^ m 
8t ' P 




dT 


(7c) 


(7d) 


Replacing the appropriate terms in equation (7a) with equations (7b) and (7d) gives 


=i_ m_9_ 

8t 8t 6 P 877 

Substituting equations (2b), (5), (6), and (8) into equation (3) gives, in the trans- 
formed moving coordinates, 

8t'\ , 1 0 ,, 8t\ 1 ® k 7^ ^ ® 


(8) 


1 


h^hs 

6^ 977 


877/ x2 9«\hl ^ HJ Xb 8^\hi ^ 6 877/ 5Xb 9T7yhi 8|y 


77A, _M^3 


8T 




= nc iI+ 5 ^£Ei 
pUt p6 877, 


(9) 
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where 


A=i-^ 

Xb 8^ 


( 10 ) 


The solution to equation (9) depends on the initial conditions and the specified 
boundary conditions. These conditions are discussed in the following sections. 


Initial Conditions 

The initial conditions that must be specified are the temperature distribution, 
mass-transfer rates, and the body shape. For most cases of interest, the initial tem- 
perature distribution is uniform and the initial mass-transfer rate is zero. 


Surface Boundary Conditions 

Two conditions must be specified at the outer surface. Either the rate of removal 
of the surface material or the surface temperature must be specified; the other condition 
is provided by an energy balance at the surface. 

S urface recession .- Ablation is assumed to result from a chemical process (oxida- 
tion) or from a phase-change process (sublimation). The rate of mass loss due to oxida- 
tion by molecular oxygen is, for an ith-order reaction. 


m 


c 


- A-c 



-Bc/Tw 

Mogj 


( 11 ) 


In this analysis, it is assumed that all oxygen at the surface is in molecular form. 

The net rate at which oxygen diffuses to the surface is, from reference 8 (assuming 
a unit Lewis number). 


“02 = 


^,net 

He - HwV ® 



( 12 ) 


where 


‘lC,net = ‘lc(l “ + ofe^s) " 



+ 



2 



(13) 
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which is the net convective heating rate to a hot ablating surface. Either transpiration 
theory (/3 = 0) or linear ablation theory (/3 = 1) can be used to account for the effects of 
mass transfer on the convective heating rate. 

The rate at which mass is removed by oxygen must be proportional to the net rate 
at which oxygen diffuses to the surface; that is, 


nif 


Xm 


02 


(14) 


Therefore, combining equations (11), (12), and (14) gives the rate of ablation by oxidation 
in a half-order reaction fi = as 


J Mw(He-Hw) 


|Mw(He - Hw) 


2 

+ 4k2c^ P ? 

2 [ ^02^C,net" | 

L '^02‘lc,net^ J 

^ ® Mq2 

J 


(15) 


where 


K = Af,e 


Bc/Tw 


The equation for a first-order oxidation reaction (i = 1) is 


m 


c - 



(16) 


Equations (15) and (16) apply to both the reaction-rate-controlled and the diffusion- 
controlled oxidation regimes as well as the intermediate conditions. 

The rate of ablation by sublimation is 


m 


s “ 




(17) 


The form of equation (17) is compatible with reference 9 and most other sublimation 
theories. Either equation (15), (16), or (17) can be the boundary condition that defines 
the rate of removal of the surface material. 
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Surface location.- The thickness at any point § at any time is 


Jq P 


( 18 ) 


The mass-loss variation about the body caused by the heating-rate distribution will 
result in a nonuniform change in the geometry of the outer surface. This effect will in 
turn affect the heating and pressure distribution as discussed in a subsequent section. 

Surface energy balance .- The heat input to the surface consists of convective and 
radiative heating and heat from combustion when oxidation occurs. This heat input must 
be accommodated by one or more of the followii^ mechanisms: (1) aerodynamic blocking 
by mass transfer, (2) reradiation from the surface, (3) conduction into the material, and 
(4) sublimation of the surface material. 


The surface energy balance is 

_ krj 9T 
^net - 5 ^ 

where 


(19) 


‘Inet = ^C,net “‘Ir + “c - mg AHg 



( 20 ) 


The heat absorbed during sublimation and the heat released during oxidation are 
considered separately in equation (20) as is the blocking effectiveness of the gases pro- 
duced by oxidation and sublimation in equation (13). In the present analysis, oxidation 
and sublimation are not allowed to occur simultaneously. 

The mass transfer affects only convective heating in this analysis. Reference 10 
indicates that at hypersonic entry velocities, radiant heating may also be significantly 
affected by mass transfer. However, at present there is no quantitative analysis for the 
effect of mass transfer on radiant heating, and it is therefore neglected. Additional 
terms can easily be included in equation (20) to accoimt for other phenomena which may 
affect the energy input to the surface. 

Correction for c hange in body geometry .- The heat input to the surface is also 
affected by changes in body geometry. In the present analysis, the stagnation convective 
and radiative heat-transfer rates are adjusted for changes in the body bluntness as 
follows: 


^ “ ^,o 



R 


stag 


and 


^r “ *lr,o 


^stag,o 

^stag 


( 21 ) 
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Geometry changes affect the heating and pressure distribution around the body as 
well. The heating-rate distribution is computed by using equation (14) of reference 11 
and modified Newtonian pressure distribution. The methods for evaluating the body 
shape parameters which are used to determine the heating-rate and pressure distribu- 
tions are given in appendix A. 

Modified NeAvtonian pressure distributions are sufficiently accurate for most appli- 
cations where the body is a hemisphere or a hemispherically blunted cone. For cones, 
the cone angle must be less than the value required to maintain supersonic flow over the 
conical portion of the body. Also, for ablating hemispherical or hemispherically capped 
bodies, the modified Newtonian pressure distribution and consequently the heating-rate 
distributions become inaccurate as the body shapes are blunted; that is, as the ratio of 
the body radius to nose radius approaches zero. 

Boundary Conditions Along the ^-Axis (0 = ^ = 1, r] = 0) 

Several options or combinations of conditions are considered along the surface 
rj = 0 (y = 0). This interior surface may be perfectly insulated or lined with a heat sink 
which may or may not radiate energy to a body at a constant temperature Tg. The 
boundary condition is 

( 22 ) 


The layer along rj - 0 is assumed to have a constant heat capacity p"Cpt". 

Boundary Conditions Along ^ = 1 , 0 = rj = 1 

The same conditions may exist along this surface that exist along rj = 0, 0 = ^ = 1. 
The boundary condition along this surface is 


^ DT 
hj 8x 


- -n'c' t' - 
” 8t 



(23) 


The thermally thin layer along 


p'c^t'. 


1 = 1 is assumed to have a constant heat capacity 


Boundary Conditions Along Line of Symmetry (4 = 0, 0 < t/ < 1) 

The assumption of an axisymmetric body requires the following boundary conditions 
along the line of symmetry: 
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(24a) 


DT _ t?A 3T 

8x 9 6 6 977 


and 


A=^-^ = 0 (24b) 

Xb94 

The boundary conditions along this line of symmetry at the surface | = 0, 77 = 1 
and at 77 = 0 are the same as those previously described. However, this line is mathe- 
matically a singularity for equation (3). The equations required along this line are dis- 
cussed in the following section. 


Singularities 

Equation (3) applies to the entire region of interest; however, a line of discontinuity 
(singularity) exists. An inspection of equation (2c) shows that at x = 0, that is, along the 
line of symmetry, the scale factor h 3 vanishes. This coordinate singularity can be 
eliminated by using proper approximations valid only near x = 0. Appendix B presents 
a detailed derivation of a form of equation (3) that applies at x = 0 and is given in 
coordinates as 


2kg/i e^T J_?7 9T^ 1 9 / 9 t\ ^ ^T 

h2 U 8|2 “ Xb 6 977 9^ 62 9t 7\J^ 97]^ 6Rhi 97] 



+ 


m77 9T 
6p dr] 


(25) 


The boundary conditions at ^ = 0 are, at 77 = 0, 


kr? 9T 
6 977 


+ erg* 

9t 


P”c”t 



(26) 


and, at 77 = 1, 


k77 9T 
6 977 


%et 


(27) 


METHOD OF SOLUTION 


The differential equations that define the temperature field in an ablating axisym- 
metric body of revolution are given in the previous section. To obtain a solution, these 
equations have been approximated by finite-difference equations and programed for solu- 
tion on a Mgh-speed digital computer. The methods used to derive the finite-difference 
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expressions are given in appendix C along with a summary of the finite-difference equa- 
tions obtained by these methods. 

The method of solution of the unknown temperature field defined by equations (C23) 
to (C31) is essentially that used in reference 12. This method is classed as an 
alternating-direction implicit method which has the advantages of being implicit, stable, 
and amenable to rapid solution. With this method, the second derivative d^T/By^ is 
replaced by a second difference evaluated in terms of unknown temperatures at time 
T + At, and the other derivative is replaced by a second difference evaluated 

in terms of known temperatures at time t. This formulation is implicit in the 
y-direction. The procedure is then repeated for a second time step of equal size, with 
the formulation implicit in the x-direction. The alternation in solution, that is, column- 
row, is continued over the specified time period. It should be noted that to maintain a 
stable numerical solution, a pair of successive row-column solutions is required. Suc- 
cessive pairs of solutions, however, may have different time steps if desired. 

Equations (C23) to (C31) take the form, for either a row or a column solution, of 

AjTj-l + BjTj . CjTj.i = Dj (28) 

where j = 1, 2, 3, . . ., S for the column solution and j = 1, 2, 3, . . ., L for the row 
solutions. Equation (28) represents L or S equations and L or S unknown tem- 
peratures. Since equation (28) results in a tridiagonal matrix of unknown temperatures, 
this set of equations can be quickly solved simultaneously by using a procedure based on 
the Gauss elimination method. This procedure is discussed in reference 12. 

RESULTS AND DISCUSSION 

In the following section, the accuracy of the solution for the unknown temperature 
field defined by equations (C23) to (C31) is evaluated by comparing the results of the 
numerical solution with two exact solutions. Also, the results obtained by application of 
the present analysis via the associated computer program to several special two- 
dimensional ablation cases are presented. Computer-drawn curves showing the com- 
puted results for some of these cases illustrate the plotting feature of the program. 

Comparison With Exact Solutions 

The exact solutions to two heat-transfer problems are used to evaluate the accu- 
racy of the numerical results. The exact solution of an orthotropic ablating body with 
temperature-dependent properties is not available. Therefore, comparisons are made 
with results from exact solutions for homogeneous, nonablating bodies with constant 
properties. 
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Insulated thick- walled hemisphere .- The exact steady-state temperature distribu- 
tion of the insulated thick-walled hemisphere shown schematically in figure 4 is derived 
in appendix D in spherical coordinates and is 


T(r,e) =:^To +-p 


Tq(s cos^ 0 - l) 


(R + 6)2 + i r5(R + 6)-3 


:2 +|-R5r“2 

O 


(29) 



The results of a comparison between the exact steady-state case and the numerical 
analysis with a 10- by 10-node network are shown in table 1. The fact that the error is 
no greater than about ±0.5 percent indicates that the errors introduced by the finite- 
difference pattern are small. 


TABLE I.- TEMPERATURE DISTRIBUTIONS OF THICK-WALLED HEMISPHERE 
FROM EXACT SOLUTION AND FROM NUMERICAL ANALYSIS 
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Percent error, — sxact ca^ ^ 100, for values of 

'T exact 

1 

O 



0 

0.125 

0.25 

0.375 

0.50 

0.625 

0.75 

0.875 

1.00 

0 

-0.162 

-0.220 

-0.234 

-0.240 

-0.243 

-0.248 

-0.257 

-0.274 

-0.298 

.25 

-.488 

-.507 

-.473 

-.407 

-.318 

-.280 

-.136 

-.088 

-.095 

.50 

-.467 

-.493 

-.463 

-.397 

-.304 

-.194 

-.093 

-.034 

-.048 

.75 

-.356 

-.395 

-.381 

-.338 

-.265 

-.168 

-.067 

-.007 

-.037 

1.00 

-.209 

-.248 

-.254 

-.245 

-.210 

-.143 

-.058 

-.003 

-.046 


One- dimensional insulated slab .- The exact transient temperature distribution in a 
one-dimensional insulated slab shown schematically in figure 5(a) is given by equa- 
tion (All) in reference 13 as (in the notation of the present paper) 
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Time, sec 


(a) Time step, O.O625 sec. 


Figure 5*- Transient temperature histories for a one-dimensional insulated slat 
= 1.13^9 Mw/m2; k = 62.35 w/m-°K; p = I60.I kg/m?; = 0.4l87 J/g-°K- 



( 30 ) 



This equation gives the temperature at a point y as a function of time. Constant ther- 
mophysical properties were used for the slab (see fig. 5(a)) and the numerical values 
were chosen to facilitate computations and not to represent any particular material. 

Comparisons between the exact solution and the numerical Solution for two differ- 
ent time steps, 0.0625 and 0.25 second, are shown in figures 5(a) and 5(b). The numer- 
ical results are in good agreement with the exact results for both time steps. 



Figure 5»- Concluded. 
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The agreement between numerical and exact results for the larger time steps is not 
as good as that obtained with the smaller time step. The maximum errors for the first 
three time steps of 0.0625 second were, respectively, 0.70, 0.36, and 0.25 percent, whereas 
the maximum errors for the larger time step were -1.22, 0.56, and 0.26 percent. In both 
cases, the solutions were stable and converged rapidly. 

Application 

The results obtained from application of the present analysis to ablating bodies are 
presented in this section. Two of the cases discussed, an ablating hemisphere-cone body 
and an ablating hemisphere, are of general interest. Two additional cases provide a com- 
parison, on a qualitative and quantitative basis, between test and analytical results. All 
the cases illustrate situations that are not amenable to a one-dimensional analysis. 

Hemisphere-cone body.- A graphite hemisphere-cone body is exposed to stagnation 
convective and radiant heating rates of 34 and 11 MW/m2, respectively, and a stream 
enthalpy of 93 MJ/kg. These energy levels are typical of earth entry at hyperbolic veloc- 
ities. The body is assumed to be subjected to this environment for 60 seconds. Fig- 
ure 6(a) shows computed geometry as a function of time during the 60-second exposure. 



(a) Shape change. 

Figure 6.- Graphite hemisphere- cone body. Initial stagnation cold-wall heating rates 
of = 34 MW/m2 and = 11 MW/mS; Hg = 93 Mj/kg; In air. 
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Note the nonuniform surface recession over the body which requires continued adjustment 
in pressure and heating-rate distributions over the body. Adjustments to the stagnation 
heating rates are reqtiired to account for increased nose bluntness. These latter adjust- 
ments, as computed with the present analysis, are shown in figure 6(b). 


stag 

or 


^'^r^'^r,o^stag 


.8 


0 


_1 1 _ 

20 4 o 


_l 

60 


Time, sec 

0>) Variation of stagnation heating rates (due to changing 
nose hluntness) with time. 

Figure 6.- Concluded. 



A 5- by 10-node network, that is, S = 5 and L = 10, was used for this example. 

A finer network would produce smoother profiles of the hemisphere-cone body than are 
shown in figure 6(a). This is particularly true in the hemispherical portion of the body. 
The total number of stations that should be used will depend on the physical size of the 
body considered and the gradients of the surface inputs, that is, heating rates and pres- 
sure on the surface. 

Hemisphere .- A low-density phenolic-nylon hemisphere, approximated with a 10- by 
10-node network, is assumed to be exposed to a convective heating rate of 3.4 MW/m2 in 
air for 110 seconds. The method presented in reference 7 was used to approximate the 
thermophysical and thermochemical properties of a charring ablator in a single-layer 
material. 

Figure 7(a) shows the outer surface location as a function of time. For these cal- 
culations the heating-rate distribution was based on Newtonian pressure distribution, and 
therefore, figure 7(a) shows qualitative rather than quantitative results for an ablating 
hemisphere. Note that for about the first 70 seconds, nonuniform ablation occurs across 
the surface; however, after 70 seconds, the recession is uniform. This behavior, although 
expected, would not be revealed in a one -dimensional analysis or a multidimensional anal- 
ysis that does not consider shape change. Similarly, the variation in the stagnation 
heating rate with nose bluntness would not be revealed in a one-dimensional analysis or a 
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(a) Shape change. 


Figure 7*- Hemispherical nose cap of lo-w-density phenolic-nylon. Initial stagna- 
tion cold-wall convective heating rate of 3.4 He = 7 Mj/kgj in air. 

multidimensional analysis that does not consider shape change. The variation in the 
stagnation heating rate with time resulting from nose bluntness, as computed by the 
present analysis^ is shown in figure 7(b). 

Teflon sphere .- The results of a series of tests on teflon spheres at various heating 
conditions are given in reference 14. One of the spheres, model E, was tested at an ini- 
tial stagnation cold-wall convective heating rate of 6.41 MW/m2, a stream enthalpy of 
4.88 MJ/kg, and a Mach number of 2.6. The experimental profile history for model E is 
shown in figure 8(a). 

The analytical profile history for model E, obtained by using the present analysis 
with a 10- by 10-node network, is shown in figure 8(b). The agreement between the 
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experimental and analytical profile histories is good, qualitatively, up to about 8 seconds. 
The analytical calculations were terminated at 8 seconds because the present analysis 
cannot consider the entire sphere, but only the forward hemisphere. Therefore, in the 
analysis the stagnation point cannot recede more than a length equal to one body radius. 
Also, the heating-rate distribution used is not accurate for highly blunted bodies. 

A comparison between the experimental and analytical stagnation-point recession 
is shown in figure 9. The analytical and experimental recessions are in good agreement. 
This agreement is better quantitatively than that shown in figure 8 because the inaccurate 
heating-rate distribution predicted with the present method for highly blunted bodies did 
not affect stagnation-point recession as directly as it affected the overall recession. 


Test data ref. 14 
Present analysis 




0 10 20 
Time, sec 

Figiire 9>- Comparison between the experimental and analytical 

stagnation-point recession for an ablating teflon sphere. 

Initial conditions: , = 6. kl MW/m?; 

(»/ ^ 

Hq = 4.88 Mj/kg; Mach number, 2.6. 

Laser test .- The results of tests of a low-density phenolic-nylon ablation material 
at high radiant heating rates produced by a laser beam are given in reference 15. The 
specimen was exposed to a laser beam, with the energy distribution shown in figure 10(a), 
for 2 seconds. The maximum radiant heating rate at the center was about 40 MW/m2. 
The final geometry of the ablated specimen was characterized by a very thin shell of 
charred material, figure 10(b). The present analysis was used to determine whether 
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sufficient radiation from the sides could limit recession along the sides and thus allow a 
thin shell of char to remain. 


The profiles computed first by assuming no radiation from the specimen sides and 
then by assuming that the sides radiate are shown in figure 10(c). Note that without radi- 
ation, the cusp profile is similar to that obtained during the test except that the sides are 
not as high. When the sides are allowed to radiate, the center recedes about the same 



(laser beam) 

i 


Estimated 
final shape 


3 .1 cm I 



cm 


(b) Estimated specimen shape. 


Side radiation 

No side radiation 



(c) Calculated recession. 

Figure 10.- Laser test results for lew-density phenolic-nylon 

specimen (ref. 15 ). 
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amount as with no radiation; however, tall sides, more characteristic of the experimental 
results, are obtained. This substantiates the conclusion drawn in reference 15 that the 
sides of the model were cooled sufficiently by radiation to permit a considerable height of 
char to remain at the edges. 

The profiles shown in figure 10(c) were computed with a 10- by 5-node network 
(i.e., S = 10 and L = 5). The computation was repeated with a 10- by 10-node network 
with essentially the same results. 

CONCLUDING REMARKS 

The differential equations governing the transient response of an ablating axisym- 
metric, orthotropic body have been derived for fixed points in a moving coordinate system. 
These equations have been expanded into finite-difference form and have been programed 
for numerical solution with an implicit technique on a digital computer. Numerical 
results compared favorably with the exact solutions for simplified conduction problems. 

The determination of the changes in the body geometry as ablation occurs, and the 
effect of these changes on the surface energy inputs, is a salient feature of this analysis. 
This feature was used to obtain satisfactory agreement between numerical and experi- 
mental results of an ablating teflon sphere and a small test specimen exposed to a high- 
intensity laser beam. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va. , April 2, 1971. 
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APPENDIX A 


SHAPE CHANGE DUE TO MASS LOSS 

Mass transfer at the surface due to ablation causes a change in material thickness. 
The change in material thickness at any point n on the surface is, from equation (18), 

6n=6n,o- (Al) 

When there are variations in heating and pressure over the body surface, the mass trans- 
fer also varies over the surface. This variation in mass transfer causes a nonuniform 
change in thickness and, hence, a change in shape which consequently alters the heating- 
rate and pressure distributions over the surface. 

The heating-rate and pressure distributions are recalculated, with the analysis of 
reference 11, as the shape changes. The methods used to evaluate the shape parameters 
""^n? "/^nj I^stag required to implement the analysis of reference 11 are given in this 

appendix. 

The w,z coordinate system shown in figure 2 is used to define the surface geometry 
at any time. The surface coordinates at station n as a function of time are 

w(t) = Rcyi + 6(t)cos 6 (A2) 


and 


z(t) = z(t)q + 6(t)q - 6(t) sin 6 (A3) 

The instantaneous value of the angle between the free-stream velocity vector 
and the local normal to the surface is required in defining a new pressure distribution. 
This angle is determined as follows: 

For n = 1, 


^1 = 0 


for 2 = n < L, 


xl/^ = 90° - tan"l 


^Wn+i - Wn-A 

^^n+l - ^n-l/ 


(A4) 


(A5) 
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and for n = L, 


= 90° - tan 


-i TL - ^L-1 ^ 
" ^L-1, 


(A6) 


Both the stagnation convective and radiative heating rates are functions of the 
instantaneous radius of curvature. The radius of curvature in the stagnation region is 
obtained by finding the radius of a circle passing through points (S,l) and (S,2) (see 
fig. 11) in w,z coordinates. 


L) 



(1,L) 


Figure 11.- Location of finite-difference stations. 

The equation of a circle in the w,z coordinate system with its center on the z-axis 
is of the form 


n 

- ^z j + Rgtag^i 


9 2 

+ w = Rgtag 


(A7) 


Evaluating equation (A7) at the point (S,2) and solving for Rgtag gives the radius of 
curvature in the stagnation region as 


R 


stag 


2 2 2 
z| - 2 z 2 Zj + z| h- 


2(z2 - zjJ 


(A8) 
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APPENDIX B 


DERIVATION OF THE DIFFERENTIAL EQUATIONS AT 
THE COORDINATE SINGULARITY 


The governing differential equation, equation (3), has a line of singularity at x = 0. 
This singularity can be eliminated by examining the behavior of the scale factor hs as 
X - 0. 


As X — 0, Rgyi — X and 


cos 6 = 


I^cyl 

R 


R’ 


Therefore, 


^3 = Rcyl + y cos e - 




= xhj 


(Bl) 


and 



,^y/x-o ^ 


(B2) 


(B3) 


Now consider equation (3) in the expanded form 

1 B /kx 8 t\ kx 8T 9^3 18/, 8 t\ % 9T ^^3 6T 

— — I + -z — T : — + — — (“iky I + 7-^ T 7 — = PCp — 


hi 8x\hi h^h ^1 

1 3 


9y/ hi 3 0y 8y 

As X — 0, equations (Bl), (B2), and (B3) reduce equation (B4) to 


-P 8 t 


(B4) 


■ ‘ ^ ^ By Rhl 8y P 8 t 


hj 0x\hi 9x j 8x hj By 


(B5) 


8T 

At X = 0, the axisymmetric-body assumption requires that - 0 and — — = 0. 

As a result, the second term on the left side of equation (B5) becomes indeterminate. 
Appl 3 dng L'Hospital's rule to this term gives 

8T 


lim 


^ 9x kx a2-p 


x-0 hjx h^ 8x2 


(B6) 
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APPENDIX B - Concluded 


Equation (B5) may now be written as 

2kx 8 / 8t\ 8T _ 8T 

hf 3x2 ay " 8 t 


(B7) 


which is the governing differential equation along the line x = 0 and is solved with the 
boundary conditions of equations (19) and (22). 
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APPENDIX C 


FINITE-DIFFERENCE EQUATIONS 


The differential equations are put in finite-difference form through the use of 
Taylor's series expansions. Forward, central, and backward differences are used. The 
methods used to obtain these differences are from reference 8. These methods, in gen- 
eral, use Taylor's series expansions at points ±A^ or ±Arj to evaluate first-order 
derivatives and Taylor's series expansions at points ±A|/2 or ±At 7/2 to evaluate 
second-order derivatives. Typical finite-difference expressions used are summarized 
in this appendix. The sketches in figure 12 illustrate the spatial relationship between 
points used in the Taylor's series expansions. 



(a) For standard derivatives. 


( Ati,x) 


2 


( AT1,x) 

+ 2 

(l,x) 


(-^,x) 1 

2 


H 

(ti, 

2 


► Ax) 


4 

(^•q,x) 

(b) For cross derivatives. 


Figure 12.- Spatial relation between points used in Taylor's series expansions. 


First-Order Derivatives 

The first-order derivatives are given by the following equations which are correct 
to order of Ar]^: 

Forward difference 


8T _ ^'^Arj ~ ^'^ri ~ '^2Ar] 

Sri 2 At? 

V 


Central difference 

8T _ ’^At7 “ '^-Af] 

dn 2 At] 

V 


(Cl) 


(C2) 
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Backward difference 

8T 3T^ - 4T_^7j + T_2A?7 
9^7 ^ ~ 2 A?7 


(C3) 


Second-Order Derivatives 

The second-order derivatives are given by the following equations: 
Forward difference 


1-B^ 
8r]\ dri 


®®An/2^'^A7j ” - ^SAr]/2(^2AT] ~ TAt?^ " ® 


3 Atj^ 


where B = hjh3k. 
Central difference 


® 1 


®Arj/2l 

~ '^rj) ~ ^-Ar]/2(^ri ~ '^-Arj) 

8rj' 


1 " 
V 

Atj^ 


(C4) 


(C5) 


Backward difference 



(C6) 


Cross Derivatives 

The methods used to evaluate cross derivatives can be shown by considering the 

term 


877 ^ 8xj 


(C7) 


Term (C7) is first expanded in the rj-direction. Thus, 


B 


-Lb^I 

8ri\ 9x^ 


9T 

9x 


- B 


Arj,x 


5T 

9x 


'-At7,x 


r],x 


2 Arj 


(C8) 
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The first derivatives are then expanded in the x-direction to give 


8 ^ 3t\ _ ®At7,x(^ At/, Ax “ "^A?7,-Ax) “ ®-At/,x('^-A77,Ax " "^-At/,-Ax) 

9r?y 4 At/ Ax 


(C9) 


Equation (C9) is correct to Ar/2, Ax^. Equations (C8) and (C9) are central-difference 
expansions. 

The forward expansion of term (C7) is obtained by evaluating at Arj and 2 Ar/, 
which gives 




4b|I 

8x 




-3b|I 

9x 


- B 


r/,x 


0T 

w 


'2An,x 


2 Ar/ 


(CIO) 


The first derivatives are then expanded to give 


8 /g 9 t\ _ ^^A7),x(^A77,Ax ~ ~ ~ '^ 77 , -Ax) ~ ^2 At;,x(^2At),Ax ~ ’^2a?;,-Ax) 

9)7 \ 9x/,j ,j 4 A )7 Ax 

The backward expansion of term (C7) is obtained by evaluating at -Ar/ and 
The resulting expression is 

8 l-rt 9t) x(Tt)|Ax “ ~ 4B-A);,x(t-A)7,Ax “ T_/^,j^_Ax) + B- 2 A) 7 ,x(T- 2 Ar 7 ,Ax “ T_2 At7,-Ax) 


(Cll) 
-2 Arj. 

(C12) 


A Special Case 

The boundaries along x = 0 and x = L require special attention. At these two 
locations the boundary conditions are, respectively. 


and 


kx DT _ /_1_ ^ TjA 8T' 
hi 9x “ ix^ 8^ " 6 8r/^ 


kx DT 
hi 8x 


_n*c' t’ 


CTc’ t: 


m,L 



(C13) 


(C14) 
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Because of symmetry at x = 0, 

A =1^ = 0 (C15) 

8x 

Therefore, at x = 0, the boundary condition becomes 


DT _ 1 ax 
9x Xjj 8^ 


(C16) 


A forward-difference expansion different from equation (C8) can now be obtained 
for the first term on the left side of equation (25). This difference expression of 82 t/9x 2 
is obtained as follows; A Taylor's series expansion of Tj^ 2 terms of Tj^^l is 


Tm,2 = Tm,l + Ax 




8x 


9x^ 


9x" 


24 


(C17) 


9x^ 


Similarly, the expansion of Tjn^S in terms of T^^ 1 is 

3 , T,,. 1 . 2 AX ^ 

, 8x ^ 


9x" 


9x" 


24 


9x" 


(C18) 


Eliminating the third derivative between equations (C17) and (C18) gives 


8T 


m,2 “ - "^Tm,! + 6 Ax 


8T 


m,l 4 Ax2 

— ^ + 


8^T 


m,l 


9x 


9x^ 


- — Ax4 
24 




0x 


(C19) 


After the boundary condition of equation (C16) is utilized, equation (C19) becomes 

^ 8Tjn,2 ~ Tm^3 - 7Tni,l 
9x^ 2 Ax2 


(C20) 


Along X = L, it is advantageous to write the first term in equation (3) in a semi- 
expanded form, which is, in the ^,77 coordinates, 






Xb9|l^ SxJ 


- T? 


A 9 4 Dt\ 

6 9?] I 8xy 


(C21) 
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when a half-station backward-difference expression (eq. (CIO)) is used, the first term on 
the right side of equation (C21) becomes 

-3B-Ax/2(-te)_^^^ 8 (C22 

Xb A{ Xb A{ 3Xb A{l,Sxjx 

The last term is evaluated from the boundary condition at x = L (eq. (C14)). The 
other terms are expressed in forward-, backward-, or central-difference form as 
required. 


Finite-Difference Approximations of the Field Equations 

The locations of the finite-difference stations (m,n) are shown in figure 11. The m 
and n subscripts correspond to the r 7 -coordinates and |- coordinates, respectively. 


The methods used to change the differential equations to finite-difference form are 
given in the preceding sections of this appendix. A summary of the finite-difference equa- 
tions obtained by these methods is presented here. 


For any point (m,n), where 


l<m<S, l<n<L, the governing equation (eq. (9)) is 





hlhs^T^ f'^m+ljn " "^m,n^ “ “ "^m-ljn^ 

o’ o’ _ 


4x^ At7 

/h3k^?7A!' 


hsk^rjAX 


h^6 


m,n+- 




n+1 “ "^m-ljn+l ■*" '^m+l,n “ "^m-ljU^ 


hjB 


^Tni+l,n “ 1’m-l,n + n-1 “ Tm-ljii-l^ 


'm,n-^ 


?7Ak^ 


V46 Ar)Xi, A|/^ 


lV^l/m+l,n' 


m+l,n+l “ "^m+ljn-lj " ( 


^m-l,n+l “ "^m-l,n-l) 


(Equation continued on next page) 
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(C23) 


At station (m = 1, n = 1), the boundary condition (eq. (26)) is combined with equa- 
tion (25) to give 

a fal,2 ~ '^1,3 ~ , 1 j^M3/2,i('^2,l ~ 




1,1 


X£ 


2 


1,1 


3 A?7 


N5/3,i('^3,1 - T2,i) _8_ 

3 A^2 3 a, 


)"c"t" + 0 E''(rf , - tJ 


At 


1,1 B 




p”c’’t'Y '^'^^’n - - t|^ 

p A At / Rhi \ 1,1 By 


(C24) 


For the curvilinear coordinate system, the a term is set equal to 2. It should be noted 
that for Cartesian coordinates, a slight modification must be made to equations (C24), 
(C25), and (C26). This modification is required since there is no line of singularity in the 
Cartesian system as there is in the curvilinear system. The correct forms of equa- 
tions (C24), (C25), and (C26) for the Cartesian system are obtained by setting the a term 
in these equations equal to 1. 

For station (1 < m < S, n = 1), equations (24) and (25) yield 



1^2. 2 
V6 Ar? ; 


N) 


m,l[ 


m+il 


m+1,1 




-Tm-l,l) 


(Equation continued on next page) 
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/ X 1 f \ 

‘ Nm,l 

For station (m = S, n = 1), equations (25) and (27) result in 

Lx^aIk^ Ws,l(®°S,2 - «S,3 - ’«S,l) 

\ b l/s,i<- ^ ^S,l 


(C25) 


+ 86 


S,1 ^’’(‘lnet)s 1 


= fpc„] ^ + — — 

\ P/S,l At 26 At7 


(sTs, 


1 + '^S-2,1 ■ '^'^S 


-l,l) 


(C26) 


For station (m = 1, 1 < n < L), equations (9) and (26) combine to give 

^h<}kt\ . /hokt 


x 2 a ^2 




2 2 
36 Atj 



3/2 


^('^2,n " '^l,n) " 


- (Shihs 


p"c”t" + o-e”(Tf - T^, 

P At \ l,n Bi 


= (hihsPCp)^ ^ 


AT 


l,n 


At 


(C27) 


For station (m = S, 1 < n < L), combining equations (9) and (20) yields: 


(hih3k^) (Ts_i,„ - V „) - 9(hih3k,) (Tg -Tg.^ ) 

G-A n G_ 1 


36^Aj7^ 


S-f,n 


S-in 


8(hih36 Ar,)g^^(q„et)^^^ 


x 2 a ^2 



S,n+1 


(^S,n+1 “ 


(Equation continued on next page) 
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(C28) 

For station (m = 1, n = L), combining equations (9), (22), and (23) gives 



(Equation continued on next page) 
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- 8 (hih 36 ^ 


1,L 


P"c^'t 


ATi T 

itfn ■*■>■*-* 


At 


+ ae" T 


1,L 



= (hlhspcp)^ ^ 


ATi^L 

At 


(C29) 


For stations (1 < m < S, n = L), equations (9) and (23) yield 



m,L-| 


3x^ A| 

+ Tm+1,L-1 “ "^m-ljE-l 
/ 77 A 


(Tm,L - Tm,L-l) ^ / tjA 













A 


(^m+l,L-l “ "^m-ljL-l "^m+l,L-2 “ '^m-l,L-2^ 


- S(h3) 


m,L 


p'c’ t' + ere' 

P At 







fr] Ak^\ 


hoTjA' 


1^1 (Tm.l,L - T„.,l) - (t„.l - Tm-1,L 

• >vi _i_ i T . ' ' _ 1 T 




m-— ,L 
2’ 




^Tm+ljL “ ~ ('^m,L “ "^m-ljE) 

TYlO-i T. w-i_ i T 


m+i,L 


m--,L 


= (hl^sPCp) 


m,L 


AT 


m,L ^ /mT 7 \ Tm+l^L ~ 


At \6p, 


m,L 


2 A?7 


(C30) 
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For station (m = S, n = L), equations (9), (20), and (23) are combined to yield 



3 A?7^6^ 


S-f,L 






S-A.L' 
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APPENDIX D 


DERIVATION OF EXACT SOLUTION 


The governing differential equation for the steady-state temperature distribution in 
an internally insulated, thick-walled, hollow hemisphere shown schematically in figure 4 
is 


1 8 


^ r2 


.2 8r\ 8r 


r^sin 9 


8 / . a 9T\ 
— sin d— — 

Be\ Be 


= 0 


(Dl) 


The boundary conditions are 


8T 

861 


61=0, 7t/2 


8T 

8r 


= 0 


(D2) 


r=R 


and 


T(R + d,e) = Tq(1 + cos26») (D3) 

The general solution to equation (Dl) by the method of separation of variables is 

(D4) 


("•-?) 

Bq Inftan + Aq 

OO 

+ ^ AnPn(cos 61) 

_C„rn + Dnr'^^-"^^ 



n=l 



where Aq, Bq, Cq, Dq, A^, 0^, and Dn are constants of integration and Pn(cos 9 ) 
is the Legendre polynomial. 

When the boundary conditions, equations (D2) and (D3), are used to evaluate the 
coefficients in equation (D4), the solution to equation (Dl) is 


T(r,0) =|To + 


Tq(3 cos 20 - l)^r2 + ^ R^r"^ 
3 j(R + 6)2 + |(R + 6)-2 r5 


(D5) 
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